% Newton-Cotes公式求数值积分
clear;
format long;
f = @(x) sin(x)./x; % 采用数组运算，使得函数可以作用在向量或矩阵上
Ck{1} = [1/2, 1/2];
Ck{2} = [1/6, 2/3, 1/6];
Ck{3} = [1/8, 3/8, 3/8, 1/8];
Ck{4} = [7/90, 16/45, 2/15, 16/45, 7/90];
Ck{5} = [19/288, 25/96, 25/144, 25/144, 25/96, 19/288];

a = 0; 
b = 1;
for n = 1 : 5 
   h = (b - a) / n;
   X = 0 : h : 1;
   Y = f(X);
   Y(1) = 1;  %  y(1)=0/sin(0)
   I(n) = sum(Ck{n} .* Y);
end
I    % 1阶到5阶Newton-Cotes求积结果
Im = quad(f, 0, 1)  % matlab求积函数